%This code carries out a single-point statisctis analysis in boundary layer
%flows using parallel computing.

%Created by Ricardo Mejia-Alvarez. Urbana, IL. 08/01/09

%modified:  01/02/2010
%           01/27/2010
%           04/06/2010
%           05/30/2010
% Julio to be more plataform independent 10/20/10


clear

% Close if there is any open pool
if matlabpool('size') > 0;
    matlabpool close force local
end
matlabpool
sizePool = matlabpool('size');

folData = input('Enter the folder in which the data is located \n','s');
if ispc==1
    folData = strcat(folData,'\');
else
    folData = strcat(folData,'/');
end

%namestart = input('enter the first part of the name of the files to process (e.g. trek00036tree.dat --> trek) \n','s');
display('1. 2D-PIV')
display('2. Stereo PIV')
kind = input('?');

if kind == 1
    namekey = '*.VEC';
else
    namekey = '*.V3D';
end

disp('Do you want to define the main direction of the flow?')
disp('That will calculate a line average and depurate your data based on residual analysis')
lineAV = input('?','s');

if lineAV == 'y' || lineAV == 'Y'
    disp('Enter the direction along which you want to line average');
    disp('1. Horizontal')
    disp('2. Vertical')
    direction = input('?');
else
    direction = 3;
end

disp('Do you want to perform your calculations over a specific region of interest? (y/n)')
ROI = input('?','s');

if ROI == 'y' || ROI == 'Y'
    display('enter the number of rows/columns you are excluding from each side (format [left right bottom top]')
    display('do this based on the way you visualize a velocity field in Tecplot')
    windowProcess = input('?');
    windowProcess = windowProcess + 1; %this set the indices at the right value
else
    windowProcess = [1 1 1 1];
end

disp('enter a correction factor for velocity scaling (1 if none exist)')
FAC = input('?');

disp('enter inversion of coordinates if necessary [x , y , z]');
disp('e.g  [-1 1 1], to invert x-axis and its associated velocity');
invertCoord = input('?');

%generates folders to store results
[folRoot,folResults,folFluc,folWz,folSw,pathJunk] = newfolder(folData,'results','fluctuations','vorticity' , 'swirlingStrength','junk');

%generates a vector with the name of the files to process
[PathFiles] = filesVector(folData,namekey);

%converts the cell array 'PathFiles' in a codistributed array and then into
%local arrays to distribute the files between the active labs.
spmd
    CodistPathFiles = codistributed(PathFiles);
    CodistPathFiles = getLocalPart(CodistPathFiles);
end

tic
display('Generating ensemble average of velocity')
spmd
    
    %ensemble average of velocity
    [Xc,Yc,~,Uxavrg,Uyavrg,Uzavrg,CHCtotal,I,J,nc,FileNames] = ensembleParallel(direction,CodistPathFiles,pathJunk,windowProcess,FAC,invertCoord);
    X = cat(3 , Xc , Yc );
    Uavrg = cat( 3 , Uxavrg , Uyavrg , Uzavrg );
    
end
[Uavrg] = Composite2GlobalEnsemble(Uavrg,sizePool); %with this line, we calculate the global ensemble average

spmd
    Uavrg = Uavrg; %#ok<ASGSL> This is done to replicate the ensemble average in all labs
end

toc

FileNames = vertcat(FileNames{:});
title = 'File names';
if ispc == 1
    saver(pathJunk,'\FileNames.dat',title,FileNames)
else
    saver(pathJunk,'/FileNames.dat',title,FileNames)
end

%pathFiles should be found again because some fields may have been deleted
[PathFiles] = filesVector(folData,namekey);

%converts the cell array 'PathFiles' in a codistributed array and then into
%local arrays to distribute the files between the active labs.
spmd
    CodistPathFiles = codistributed(PathFiles);
    CodistPathFiles = getLocalPart(CodistPathFiles);
end

Il = I{1};
Jl = J{1};
ti = num2str(Il);
tj = num2str(Jl);

tic
display('Generating fluctuation fields and Reynolds stresses')
spmd
    
    %this function calculates the ensemble average of Reynolds stresses,
    %Swirling Strength, and Vorticity. In addition, it calculates and saves
    %instantaneous fields of fluctuations, vorticity, and swirling strength.
    [RestressAv,WzAv,WzRMS,SwAv]...
        = ReStressVortSwirl(CodistPathFiles,folFluc,folWz,folSw,Uavrg,windowProcess,FAC,invertCoord);
    WzRMS = WzRMS.^2;
    
end

%finding paths for swirling strength fields
if ispc == 1
    pathSw = strcat(folSw , '\');
else
    pathSw = strcat(folSw , '/');
end
[SwPathFiles] = filesVector(pathSw,'*.dat');

%converts the cell array 'SwPathFiles' in a codistributed array and then into
%local arrays to distribute the files between the active labs.
spmd
    SwCodistPathFiles = codistributed(SwPathFiles);
    SwCodistPathFiles = getLocalPart(SwCodistPathFiles);
end


spmd
    %this function calculates the RMS of swirling strength
    SwAvaug = cat(3 , zeros( J , I ) , zeros( J , I ) , SwAv);
    [SwRMS] = RMScalculator(SwCodistPathFiles,SwAvaug);
    SwAvaug = [];
    SwRMS = SwRMS( : , : , end);
    SwRMS = SwRMS.^2;
end

toc

tic
display('Saving results')

%obtaining global ensembles of all single point statistics from their
%composite arrays. This is conserving the matlab matrix format

[X] = X{1};
Uavrg = Uavrg{1};   %this is a replicated composite array
[RestressAv] = Composite2GlobalEnsemble(RestressAv,sizePool);
WzAv = Composite2GlobalEnsemble(WzAv,sizePool);
SwAv = Composite2GlobalEnsemble(SwAv,sizePool);
WzRMS = Composite2GlobalEnsemble(WzRMS,sizePool);
WzRMS = sqrt(WzRMS);
SwRMS = Composite2GlobalEnsemble(SwRMS,sizePool);
SwRMS = sqrt(SwRMS);

%Creating a matrix that defines the valid points in the ensembles.
[CHCtotal] = Composite2GlobalSum(CHCtotal,sizePool);
[CHCtotal] = depurationCHC(CHCtotal);

%Obtainig wall-normal profiles of single point statistics
if direction ~= 3
    
    if direction == 1
        Yc = X(: , 1 , 2);
        
        CHCtotalc = sum(CHCtotal,2);
        SwAvc = sum(SwAv,2) ./ CHCtotalc;
        WzAvc = sum(WzAv,2) ./ CHCtotalc;
        SwRMSc = sum(SwRMS,2) ./ CHCtotalc;
        WzRMSc = sum(WzRMS,2) ./ CHCtotalc;
        
        CHCtotalc = repmat(CHCtotalc , [1 , 1 , 3]);
        Uavrgc = sum(Uavrg,2) ./ CHCtotalc;
        
        CHCtotalc = repmat(CHCtotalc , [1 , 1 , 2]);
        Restressc = sum(RestressAv,2) ./ CHCtotalc;
        
        curves = [Yc(: , 1 , 1) , Uavrgc(: , 1 , 1) , Restressc(: , 1 , 1) ,...
            Restressc(: , 1 , 2) , Restressc(: , 1 , 3) , Restressc(: , 1 , 4) ,...
            Restressc(: , 1 , 5) , Restressc(: , 1 , 6) , WzAvc , WzRMSc ,...
            SwAvc , SwRMSc];
        
    else
        Yc = X(1 , : , 1);
        Yc = Yc';
        
        CHCtotalc = sum(CHCtotal,1);
        SwAvc = sum(SwAv,1) ./ CHCtotalc;
        WzAvc = sum(WzAv,1) ./ CHCtotalc;
        SwRMSc = sum(SwRMS,1) ./ CHCtotalc;
        WzRMSc = sum(WzRMS,1) ./ CHCtotalc;
        
        CHCtotalc = repmat(CHCtotalc , [1 , 1 , 3]);
        Uavrgc = ( sum(Uavrg,1) ./ CHCtotalc )';
        
        CHCtotalc = repmat(CHCtotalc , [1 , 1 , 2]);
        Restressc = ( sum(Restress,1) ./ CHCtotalc )';
        
        curves = [Yc(1 , : , 2) , Uavrgc(1 , 1 , 2) , Restressc(: , 1 , 2) ,...
            Restressc(: , 1 , 1) , Restressc(: , 1 , 3) , Restressc(: , 1 , 4) ,...
            Restressc(: , 1 , 6) , Restressc(: , 1 , 5) , WzAvc , WzRMSc ,...
            SwAvc , SwRMSc];
    end
    
    
    title = strcat('VARIABLES= "y", "U", "uu", "vv", "ww", "uv", "uw", "vw", "Wz", "WzRMS", "Swirl", "SwirlRMS"');
    if ispc == 1
        name = strcat(folResults,'\profiles.dat');
    else
        name = strcat(folResults,'/profiles.dat');
    end
    dlmwrite(name,title,'delimiter','','newline','pc')
    dlmwrite(name,curves,'-append','newline','pc');
    
end

%saving contour maps
Uavrg = cat(3 , X , Uavrg);
[Uavrg] = mixingMat(Uavrg);

RestressAv = cat(3 , X , RestressAv);
[RestressAv] = mixingMat(RestressAv);

WzAv = cat(3 , X , WzAv);
[WzAv] = mixingMat(WzAv);

SwAv = cat(3 , X , SwAv);
[SwAv] = mixingMat(SwAv);

WzRMS = cat(3 , X , WzRMS);
[WzRMS] = mixingMat(WzRMS);

SwRMS = cat(3 , X , SwRMS);
[SwRMS] = mixingMat(SwRMS);

if ispc ==1
    slash = '\';
else
    slash = '/';
end

title = strcat('VARIABLES= "X", "Y", "U", "V", "W", ZONE',' I=',ti,' J=',tj,' F=POINT');
saver(folResults,[slash 'Uavrg.dat'],title,Uavrg)

title = strcat('VARIABLES="X", "Y", "uu", "vv", "ww" ,"uv", "uw", "vw"',' ZONE',' I=',ti,' J=',tj,' F=POINT');
saver(folResults,[slash 'ReStress.dat'],title,RestressAv)

title = strcat('VARIABLES= "X", "Y", "<greek>W</greek><sub>Z</sub>", ZONE', ' I=',ti,' J=',tj,' F=POINT');
saver(folResults,[slash 'VorticityZ_Avrg.dat'],title,WzAv)

title = strcat('VARIABLES= "X", "Y", "<greek>l</greek><sub>ci</sub>", ZONE', ' I=',ti,' J=',tj,' F=POINT');
saver(folResults,[slash 'SwirlingStrength_Avrg.dat'],title,SwAv)

title = strcat('VARIABLES= "X", "Y", "<greek>W</greek><sub>Z,RMS</sub>", ZONE', ' I=',ti,' J=',tj,' F=POINT');
saver(folResults,[slash 'VorticityZ_RMS.dat'],title,WzRMS)

title = strcat('VARIABLES= "X", "Y", "<greek>l</greek><sub>ci,RMS</sub>", ZONE', ' I=',ti,' J=',tj,' F=POINT');
saver(folResults,[slash 'SwirlingStrength_RMS.dat'],title,SwRMS)

toc

matlabpool close



